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ABSTRACT 

Context. New high-resolution spectropolarimetric observations of solar prominences require improved radiative modelling capabili- 
ties in order to take into account both multi-dimensional - at least 2D - geometry and complex atomic models. 
Aims. This makes necessary the use of very fast numerical schemes for the resolution of 2D non-LTE radiative transfer problems 
considering freestanding and illuminated slabs. 

Methods. The implementation of Gauss-Seidel and successive over-relaxation iterative schemes in 2D, together with a multi-grid 
algorithm, is thoroughly described in the frame of the short characteristics method for the computation of the formal solution of the 
radiative transfer equation in cartesian geometry. 

Results. We propose a new test for multidimensional radiative transfer codes and we also provide original benchmark results for 
simple 2D multilevel atom cases which should be helpful for the further development of such radiative transfer codes, in general. 
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Our primary aim is thus to improve diagnostics based, in 
particular, on He i lines by treating a detailed He-atomic model 
including the atomic fine structure together with 2D non-LTE 
radiative transfer. We are motivated here by new solar promi- 
nences observations (see e.g., Paletou et al. 120011 Merenda et 
al. 120061) which also triggered some revisions of inverting tools 
(Lopez Ariste & Casini 120021) . And up to now, the later diag- 
nostic tools are limited to the assumption that the relevant (ob- 
served) He spectral lines are optically thin; it is, however, easy to 
check from high spectral resolution observations that a spectral 
line like D3 of He 1 in the visible, for instance, is not always opti- 
cally thin, even in quiescent prominences (Landi Degl'Innocenti 
1982, Lopez Ariste & Casini 2002). Furthermore, the expected 
optical thicknesses of 1 to 10 say, in such structures let us fore- 
cast the presence of significant geometrical effects on the mech- 
anism of formation of this spectral line, as the ones already put 
in evidence on Ha by Paletou (119971) . 

The combination of 2D geometry with a very detailed atomic 
model for He obviously requires a more efficient radiative trans- 
fer code as compared to the one developed from the MALI 
method by Paletou (1995). And clearly enough, the planned 
improvement of the radiative modelling including multidimen- 
sional geometry together with multi-level, realistic atomic mod- 
els have to rely on those new and fast radiative transfer methods 
based on GS/SOR with multigrid numerical schemes. 

GS/SOR methods, best implemented within a short charac- 
teristics formal solver, have been described in every details by 
Trujillo Bueno & Fabiani Bendicho (1995 ) but only in the frame 
of the two-level atom case and in ID geometry. In another ar- 
ticle, Fabiani Bendicho et al. (1997 ) have nicely described the 



1. Introduction 

Efficient iterative schemes have been introduced in the field of 
two-dimensional (2D) non-LTE numerical radiative transfer dur- 
ing the last, say, fifteen years. These developments most of- 
ten rely on the combination of the short characteristics (SC) 
method for the so-called formal solution of the radiative trans- 
fer equation (in cartesian geometry see e.g., Kunasz & Auer 
1988] Auer & Paletou 119941 and in various other geometries, 
van Noort et al. 120021) and efficient iterative schemes such as 
Gauss-Seidel and successive over-relaxation (GS/SOR) iterative 
processes (Trujillo Bueno & Fabiani Bendicho 1995, Paletou & 
Leger l2007l) together with multi-grid (MG) methods (Auer et al. 
[19941 Fabiani Bendicho et aL [19971 

Hereafter, we are interested in a more realistic modelling of 
isolated and illuminated structure, such as prominences hang- 
ing in the solar corona (see e.g., Paletou 1995, 119961) . Our fu- 
ture work will emphasis on the synthesis of the H and He spec- 
tra. Indeed, the most recent modelling efforts concerning the 
synthesis of the H and He spectrum in prominences has been 
performed using either mono-dimensional (ID) slabs (Labrosse 
& Gouttebroze 120011 [2004), or using 2D cartesian slabs in 
magnetohydrostatic equilibrium and the Multilevel Accelerated 
Lambda Iteration (MALI) technique for the solution of the 
non-LTE transfer problem (Heinzel & Anzer [200T1 I2QQ5T) . 
Concerning the spectrum of H, Gouttebroze (2006) also pre- 
sented promising new results using 2D cylindrical models of 
coronal loops. 
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Fig. 1. Schematic view of the short characteristics method for 
the computation of J in a ID grid of N points in the case of a 
ALI iterative scheme. During the first (downward here) pass the 
angular integration as in Eq. ^ is made over half an hemisphere; 
it is completed during the 2nd (upward) pass, allowing for the 
source function update according to Eq. ([T]). 



implementation of non-linear multi-grid techniques, using an ef- 
ficient iterative method such as a GS/SOR scheme. However, 
they just did not describe the implementation of, ID or multi- 
dimensional, multilevel GS/SOR scheme using the SC method. 
Besides, Paletou & Leger (120071) have finally made explicit the 
implementation of GS/SOR iterative schemes in the multi-level 
atom case, restricted though to a ID plane-parallel geometry. 

The present article aims therefore at "filling the gap" by pro- 
viding all the elements required for a successful implementa- 
tion of a GS/SOR iterative scheme in a 2D cartesian geometry. 
In order to do so, we adopt the line of detailing the method in 
the frame of the 2-level atom given that our detailed descrip- 
tion of the multilevel strategy published elsewhere (Paletou & 
Leger [20071) does not need to be commented any further for the 
jump from ID to 2D. Therefore, we also provide hereafter vari- 
ous benchmark results for the 2D -multilevel atom case, still un- 
published to date, using simple atomic models taken from Avrett 
(1968; see also Paletou & Leger 2007). Moreover, an original 
comparison between 2D numerical results and independent an- 
alytical solutions is made. 

We shall recall in §2 the basic principles of ALI and GS/SOR 
iterative schemes in the frame of a two-level atom model and 
in ID geometry. Then in §3 we shall describe, in details, how 
the GS/SOR numerical method can be implemented for the 
case of 2D slabs in cartesian geometry, therefore upgrading the 
2D short characteristic method initially published by Auer & 
Paletou (119941) . A new test for numerical radiation transfer codes 
is briefly presented in §4. Then we shall finally present, in §5, 
benchmark results for simplified multilevel atomic models in 2D 
geometry and some illustrative examples clearly demonstrating 
in which conditions geometrical effects should be seriously con- 
sidered. 



2. Gauss-Seidel and SOR iterative schemes basics 

In the two-level atom case, the non-LTE line source function, 
assuming complete redistribution in frequency, is usually written 

as 



where r is the optical depth, S * is the thermal source function 
and s is the collisional destruction probability; unless explicitly 
mentioned, S * = sB, where B is the Planck function. J is the 
usual mean intensity defined as 



J 4tt 
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(2) 



where the optical depth dependence has been omitted for the 
sake of simplicity; as usual, I v q is the specific intensity and cp v is 
the line absorption profile. Usually again, the mean intensity is 
written as the formal solution of the radiative transfer equation 
i.e., 



/= A[S]. 



(3) 



Following the Jacobi-type it erative scheme introduced in nu- 
merical transfer by Olson et al. (1986 ), we shall consider a split- 
ting operator A* equal to the exact diagonal of the true operator 
A. Now introducing the perturbations 



A = A* + (A - A*) 

5 (new) = 5 (old) + SS 



(4) 



in Eq. ([T]), we are led to an iterative scheme such that 

^(new) = [i _ (i _ s)A*] _1 {(l - e)(A - A*)S (old) + sB}. (5) 

Running the later scheme to convergence is better known in nu- 
merical radiative transfer as the 'ALI method". As schematized 
in Fig. [T] using the short characteristics method in ID geom- 
etry (Olson & Kunasz [T9871 Kunasz & Auer [T988l the for- 
mal solution is obtained by sweeping the grid say, first in di- 
rections -Q (p < 0) i.e., from the surface down to the bottom 
of the atmosphere, and then in the opposite, upward directions 
+Q (p > 0) starting from the bottom of the atmosphere up to its 
surface though. The specific intensity 7 v n is then advanced step 
by step during each pass, partially integrated over angles and 
frequencies during the downward pass while, during the second 
(upward) pass, completion of the angular integration allows for 
the full determination of the mean intensity Jk at each depth i> 
of the ID grid and, finally, to the update of the source function 



-» (new) 



S[ old) + AS k , 



k k 

on the basis of increments such that 



AS k = 



(l-s)j£ ld) +sB k -Sf d) 

l-(l-£)Ajfc* 



(6) 



(7) 



where A^ is a scalar equal to the diagonal element of the full 
operator A at such a depth in the atmosphere and where super- 
scripts (old) denote quantities already known from the previous 
iterative stage. 

For a Gauss-Seidel iterative scheme, the sweeping of the at- 
mosphere is identical but as soon as the mean intensity Jk is fully 
computed at depth-point k in the atmosphere during the upward 
pass, the local source function is updated immediately i.e., be- 
fore completion of the 2nd pass, using increments which have 
now turned into 



AS 



(GS) _ 



ax 7(old and new) , n 
- £)J) +sB k 



, (old) 



1-(1-*)A ( 



(8) 



kk 



S(T) = (l-fi)/(T)+S*(T), 



where the quantity / 



7(old and new) 



means that at the spatial point 



(1) k the mean intensity has to be calculated via a formal solution 
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Fig. 2. The 2D grid is swept four times: first pass in directions \ 
defined in panel 1 , second pass in directions \ defined in panel 
2, third pass in directions /* defined in panel 3 and fourth pass 
in directions / defined in panel 4. Sweeping must be done away 
from the boundaries so that upwind intensities - see Eq. ( [TT] ) - 
are always known. 



of the transfer equation using the "new" source function values 
^(new) already obtained at points j = ND, ...,(£+ 1) and the 

"old" source function values 5^ old) at points j = k, (k - 1), . . . , 1. 

Finally, as a next step SOR iterations can simply be imple- 
mented using 



AS 



(SOR) 



qjAS 



(GS) 



(9) 



where co is an over relaxation parameter such that 1 < co < 2. 
For two-level atom models in ID, this method was o riginally 
proposed by Trujillo Bueno & Fabiani Bendicho (119951 . 



3. The 2D-cartesian geometry case 

Hereafter, we shall describe in every details how the GS/SOR 
numerical method can be implemented for the case of 2D free- 
standing slabs modeled in cartesian geometry. 



3.1. SC in 2D: an overview 

We shall initially follow and therefore upgrade the formal solver 
of reference proposed originally by Kunasz & Auer (1988 ) and 
modified by Auer & Paletou (119941) . 

Using SC in 2D geometry, the formal solution is obtained by 
sweeping the grid four times, as schematized in Fig.[2j say first 
increasing y and z i.e., along directions £l\ (note that z = is the 
surface of the atmosphere), second decreasing y and z along di- 
rections £2 2 , third increasing y and decreasing z along directions 
£^3, and finally decreasing y and increasing z along directions 
Q4. The specific intensity I v q is therefore advanced step by step 
during each pass, partially integrated over angles, quadrant af- 
ter quadrant, and over frequencies during the first three passes 
while, during the fourth pass, the mean intensity / can be fully 
computed, completing therefore the numerical evaluation of the 
formal solution 



J(ij) - A(/j)[S] . 



(10) 



Except at the boundary surfaces where the incident radiation 
is known a priori, along each direction the specific intensity at 
the inner grid points is advanced depth after depth. As displayed 
in Fig. [3] the short characteristic starts at grid point o (/, j) and 
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Fig. 3. Example of a short characteristic across a 2D cartesian 
grid at depth point o (i, j) for a ray propagating along directions 
\ defined in panel 1 of Fig. [2j three points a, b and c are used 
for monotonic parabolic interpolation in z, in order to evaluate 
quantities at point u following Auer & Paletou ( 19941) . 



extend in the "upwind" and "downwind" directions until it hits 
one of the cell boundaries either at point u or at point d that 
is, not grid points in general. The specific intensity is therefore 
computed, according to Kunasz & Auer (1988), as 



I = I u e~^ + ^ U S U + V S + V d S d . 



(11) 



where the first part of the right-hand side of this expression cor- 
responds to the part transmitted from the "upwind" point u down 
to the current point 0, and the three last terms result from the an- 
alytic integration of 



I 



S(r)e- T dr 

Jo 



(12) 



along the short characteristic going from u to o\ expressions for 
the ^Fs can be found in Paletou & Leger (2007]). 

As shown in Fig. |3j in 2D geometry, I u , S u and Sj are not 
grid points, and they must be evaluated by interpolation on the 
basis of a set of grid point. In order to do so, one has first to de- 
termine on which axis, y or z, the upwind and downwind points 
shall lie. We introduce c y (respectively c z ) the cosine between 
the direction into which the photon is moving and the j-axis (re- 
spectively the z-axis), Ay the y length of the cell containing both 
u and o grid points, and Az its length in z. If 

Ay Az 

Cy C z 

the ray hits the j-axis and Ar u = Az/c z . 

Following Auer & Paletou (1994 ), I u and S u are determined 
by interpolation along the upwind grid-line passing through 
points a and b. To perform a parabolic interpolation, we shall 
therefore use three grid points a, b and c as displayed in Fig. [3] 
and where "quantities" have already been updated; along z-lines, 
interpolation weights would be given by 
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Fig. 4. Let i be the index along the y-axis and j along the z-axis; 
we consider the specific intensity evaluation at any inner grid 
point (ij) during the fourth pass in the 2D grid, corresponding 
to the directions D4 defined in Fig. [2] At this stage, all "new" 
grid points have already been swept, so that source functions at 
these points have been updated too. 



and similar weights should be used for interpolation in y, using 
grid points (i - (i - IJ - 1) and (i - IJ - 2) though. Then, 
we are able to calculate the upwind specific intensity as 



ceding steps. Indeed, using an expression similar to the one in 
Eq. ([T4|), for an interpolation along the z-axis, we would have 



c (new) y ^ (new; y ^ 

a,(ijr 0+1 J) ^ ^bXijf (i+hj-D ^ ^cXijr (i+U-2) 



i (new) 



1 (new) 



. (16) 



Before integrating over all frequencies and over the angles corre- 
sponding to the Q4 directions in order to obtain the partial mean 
intensity 



r dn r 



(f)yIyQdV , 



(17) 



we shall have to correct the specific intensity calculated during 
the first three passes for consistency with the source function 
updates. More specifically, the term was calculated during 
the third pass a^] 



/(old) 



/ c /'(OLD) 



"0 



^d S d 



(18) 



using 5'^ (0LD) instead of the new value s^ (new) obtained from 



the interpolation using the updated points (/ + IJ), (/ + IJ - 1) 
and (i + 1 J - 2) as shown in Fig. [4] since we have the identity 
; a correcting term 



c/ _ o/ 
d ~ u 



Aj £j) = £ f K* ew) - <H d9) 



(14) 



where specific intensity values have already been computed at 
grid points a, h and c. This is guaranteed by sweeping the grid 
away from one of the upwind boundaries. Note also that S u and 
Sd are also evaluated from (S a , St, S c ) using similar expres- 
sions. 

For the sake of accuracy and in order to avoid the genera- 
tion of spurious upwind intensities by high-order interpolation, 
one must use a monotonic interpolation i.e., set I u (and Id) equal 
to the minimum or maximum of I a and 1^ if the parabolic inter- 
polant lies outside the interval [min(4,4), max(/ fl ,/&)], as pro- 
posed by Auer & Paletou (119941) . 



must therefore be added to the the total mean intensity by in- 
tegrating the specific intensity correction over frequencies and 
over all angles Q3 - see Fig.|2j this step is equivalent to the com- 
putation of the AJ™ correction mentioned by Trujillo Bueno & 
Fabiani Bendicho fl995t in their Eq. (39). 



The two other terms /X and L\ calculated during the first 
and the second passes are also still inconsistent with the last 
source function updates because they were calculated as 



A 



(i,J) 



and 



r\(OLD) _Ar> , vtAc\(OLD) 

1 U e T It/ O h 



+ *p\s Vold) + *f\s V0LD) 

d d 



(20) 



3.2. Implementation of GS/SOR In 2D 

Assume that one has already swept the grid three times as de- 
scribed in Fig. [2] By analogy with the GS/SOR numerical strat- 
egy in ID geometry, we are now going to update the source 
function at each grid point during the fourth pass of the SC- 
2D scheme, according to the correction given in Eq. ^ and 
before passing to the next depth point. It is a quite straightfor- 
ward task at the boundary surfaces since the incident radiation 
field is known a priori from the (given) external conditions of 
illumination. 

We shall hereafter describe what has to be done at the inner 
grid points. Fig. [4] describes the situation once arriving at (ij) 
after the 2D grid was swept thrice. Using superscripts defined in 
Fig.[2j the current specific intensity comes from 



T u _|_ vp/ c /(new) + vp/ c /(old) + vp/ e /(old) 
u u d d 1 



(15) 



where one must understand quantities with superscripts (new) 
such as resulting from interpolations along upwind grid lines, 
using source functions that has been obtained during the pre- 



'/ 



/ \(OLD)^_Ar i > + *|A s \(OLD) 



+ y'Ns Vold) + yNs V0LD) 

d d 



(21) 



where we have the following identities S^ (0LD) = £^ (0LD) and 



1 \(OLD) 



, \(OLD) 



These (OLD) source functions could now be calculated using 
updated values. For example the new value sA (new) = £\( new ) j s 
obtained from an equation similar to Eq. (14) with an interpola 
tion along j-axi^j using S^f^ - and one can see, using Fig. |j 
that (ij - 1) is a "new" grid point whereas (i - IJ - 1) an 
(i - 2J - 1) are "old" grid points i.e., 



o \(new) _ \ 

J « ~ «,(ur (u-i) 



_ v ^(new) \ ^(old) \ ^(old) n0 x 



1 We use the superscript (OLD) to emphasize terms which need to be 
replaced by new values according to updated points in Fig. [4] whereas 
(old) terms remain unchanged. 

2 For an interpolation along z-axis, there are no "new" grid points to 
consider. 
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Table 1. Computation time (for a Pentium-4 @ 3 GHz processor) and number of iterations for the H i multilevel benchmark model 
of Avrett (1968) in a 2D grid with y max = 5 000 km and z m ax = 30 000 km together with 3 angles per octant and 8 frequencies; the 
temperature of the atmosphere is 5 000 K and the gas pressure p g = 1 dyn cm -2 . 



Points number 


MALI 2D 


GSM 2D 


SOR 2D 


MG 2D 


Rc 


123x123 
163x163 
203x203 
243x243 


3min9s (46) 
9min39s (79) 
22min47s (116) 
45min32s (158) 


2minl9s (29) 
6min56s (48) 
14min36s (68) 
29minl0s (90) 


lminl7s (16) 
3min33s (24) 
7min34s (33) 
14min3s (43) 


55s (11) 
lmin52s (13) 
2min50s (14) 
4minl3s (14) 


1.1 x 10- 2 
2.1 x 10" 3 
5.7 x 10- 4 
1.9 x 10- 4 



Similarly, the new value £^ (new) = ^\( new ) j s obtained using 



an interpolation along y-axis, for instance, involving S and 

S [+2^+1 ~ wim m is time, using Fig. 51 grid points at (i + 1 J + 1) 
and (i + 2, j + 1) are "new" whereas + 1) is an "old" grid point 
i.e., 



£\(new) _ \ ^r(old) 



_i_ , .\ c (new) 



+ (23) 

c,(i,j) 0+2 v 7 



By analogy, old specific intensities 7^ (0LD) and 7^ (0LD) must 
be updated to obtain new values calculated with interpolations 
using "new" grid points. 



We shall then have to calculate two other corrections AJ, 

A 



\ 
Oj) 



and AJ^ by integrating these corrected specific intensities over 
frequencies and over directions Di and Q2, following an equa- 
tion similar to Eq. (19). Finally we shall add three correcting 
terms to compute the correct total mean intensity at the current 
grid point 



J (i,j) ^ J (i,j) ^ J (i,j) ^ J (i,j) 



+ AjS + ax\ + A/\ 

Oj) 



(24) 



Then it is straightforward to update the local source function 



4 n / w) viaEq.(|8j 

However, before advancing to the next depth point (ij +1), 
it is important to add the following corrections to the specific 
intensities of the three first passes, due to the source function 
update which has just been made at the current depth point : 



( A/ A 

0j) 



0j) 



w\ [0(111 
T o Po, 



(new) _ p(old)l 

%i) a J) J 



(new) 
j) 



, (old) 



(25) 



0j) L 0j) (ij) 1 



This last stage is analogous to th e corr ection described by 
Trujillo Bueno & Fabiani Bendicho (1995 ) in their Eq. (40). 

Finally, a two-dimensional SOR iterative scheme is built 
when, at each depth-point (i, j), the source function is updated 
according to 



AC (SOR) _ , lA c(GS) 



(26) 



where oj is computed exactly in the same way as in the ID case. 



3.3. Additional notes on the whole numerical scheme 

As in the ID case, implementing a GS/SOR solver requires to 
properly order the various loops; starting from outer to inner 
loop one may find: (1) the directions Q; as shown on Fig. [2] (2) 



the direction cosines in each quadrant Q/ and, finally (3) the fre 
quencies. The corrections described in Eqs. ( [T8] ), p0| ) and (21 



require some bookkeeping of variables such as all the ^ M 's and 
the ^/s computed during the three first passes (for the further 
computation of the mean intensity). 

Details upon the implementation of GS/SOR for multilevel 
atom models were given by Paletou & Leger ('2007'). The main 
difference with the two-level atom case is the propagation of the 
effects of the local population update: it generates for each al- 
lowed transition changes in the absorption coefficients at line 
center and in the line source functions. 

Furthermore, we have also embedded the above-described 
2D-GS/SOR scheme into a nested multigrid radiative trans- 
fer method following the precise description given by Fabiani 
Bendicho et al. (1997 ). We use three grids with a grid-doubling 
strategy. On the coarsest grid (i.e., level / = 1), we iterate to 
convergence i.e., until R c i.e., the relative error on the level- 
populations from an iteration to another is "small" using the 2D- 
GS/SOR scheme. For each grid / = 2, 3 where grid level / = 3 is 
the finest one, we interpolate populations onto grid level / using 
those obtained onto grid level (/ - 1) and calculate the corre- 
sponding absorption coefficients and source functions. We iter- 
ate onto grid level / using the standard multigrid method from 
grid level / down to grid level 1=1 only until the following 
stopping criterion is satisfied 



R c (iter, I) 



I- A 



< ^R c (iter =1,1) 



(27) 



where X ^ R c (iter, l)/R c (iter - 1,1), as proposed by Auer et al. 
(TT9941) . 

We remind here the main steps of one standard multigrid it- 
eration: make one pre-smoothing iteration onto grid level / using 
a pure GS iterative scheme, then a restriction down to grid level 
/ = 1 to compute the coarse-grid equation, solve the coarse-grid 
equation onto grid level / = 1 using the 2D-SOR scheme, make 
a prolongation up to grid level / to obtain a new estimate of the 
populations, then one post-smoothing iteration onto grid level / 
using again a pure GS iterative scheme (it is important to note 
that one must make one pre- or post-smoothing iteration on each 
grid level using a pure GS iterative scheme). We used a cubic- 
centered interpolation for the prolongation and the adjoint of a 
nine -point prolongation for the restriction (see e.g., Hackbusch 
1985). 



4. Validation vs. an analytical solution 

There is no analytical solution for 2D non-LTE radiative trans- 
fer. However it is possible to compare 2D numerical solutions to 
ID solutions for which accurate and robust numerical and ana- 
lytical methods exist. In order for this comparison to be accurate, 
the slab has to be sufficiently extended in the y direction i.e. "ef- 
fectively" infinite. 
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Fig. 5. Relative errors between spatial averages of the angular 
moments /, H and K given by the 2D-SOR 2-level iterative pro- 
cess and their analytical values for e = 0.01 vs. the number of 
spatial points of a square 2D grid of extension r* = 100 in each 
direction. 



We have used the ARTY code for the computation of refer- 
ence, analytical solutions (Chevallier & Rutily 2005; see also the 
Appendix) obtained using the method of the finite Laplace trans- 
form. This code can solve, indeed, standard ID problems with an 
intrinsic accuracy better than 10" 10 ; it has already been useful in 
order to test the ALI method plus a SC formal solver in ID for 
the case of a non-illuminated, homogeneous and isotropic plane- 
parallel slab with internal, homogeneous sources (Chevallier et 
al. 2003). 

A stringent test for our 2D code was to consider a point- 
source located at the center of a non-externally illuminated and 
homogeneous slab. The central source emits isotropically in 
space, which in fact corresponds to a line, infinite along x, of 
sources. This idealized model captures most of the difficulties 
met by the numerical methods to solve the radiative transfer 
equation: the scattering is not neglected - it can also be dominant 
-, and the ponctual source will lead to large gradients difficult to 
handle when dealing with discretization of the slab. 

Then we have computed the properties of the radiation field 
emerging at the top surface of the slab (z = 0) at one frequency, 
as described by the usual first three moments /, H , K of the spe- 
cific intensity; more precisely, the later were integrated in space, 
along the top surface of the slab, in order to be compared to the 
ID analytical solutions. 

To achieve this test, we have chosen the difficult case of a 
slab of optical thickness r* = 100 in both directions where scat- 
tering dominates the absorption adopting the value s = 0.01. 
This medium is therefore effectively thick because the thermal- 
ization depth I « 1 / V3e is much less than the optical thickness 
in such as case. We used Carlson's "Set A" (1963 ) with 10 points 
per octant to describe the angular dependence of the radiation 
field and only one frequency-point. The Dirac thermal source 
term has been modelized by a sharp, normalised 2D-Gaussian 
function having half-width at half-maximum 0.16 in r. The grid 
is logarithmically refined near the center of the slab in order to 
describe accurately the shape of the 2D-Gaussian: the closest 
points from the center are at a distance 10" 4 along the axes, in 
order to accurately describe the gaussian shape whose numerical 
integral over the space has to be the closest as possible to unity. 

The two-level 2D-SOR iterative process was iterated until 
convergence of /, H and K i.e., when the second digit of their 



Fig. 6. Rates of convergence for the SOR (solid lines) and SOR- 
MG (dashed lines) 2D multilevel iterative processes. A spatial 
grid of 163 points per direction with z ma x - 30000 km and 
y m ax = 5 000 km; the temperature of the atmosphere is 5 000 
K and the gas pressure fixed at p g = 1 dyncm -2 . For bench- 
mark purposes, we adopted the simplified 3 -level H i atom model 
taken from Avrett (119681) . 



relative error did not show any more variation from one itera- 
tion to the other. For such a case, 500 iterations are sufficient. 
In Fig. [5] we demonstrate how these errors behave with the re- 
finement of the spatial quadratures; the absolute values of the 
reference solutions are given in Appendix, as well as the source 
functions and values of the specific intensity in the directions 
corresponding to the angular quadrature chosen here. 

The important point to rise here concerning this new test is 
that (i) acceptable relative errors, say better than 5% are obtained 
only for very refined grid which (ii) can hardly be handled us- 
ing a simple Jacobi-like iterative scheme such as ALI. This jus- 
tify again the adoption of very high rate of convergence methods 
such as GS/SOR plus MG. Finally, we are conducting more com- 
prehensive tests of this nature which results will be published 
elsewhere. 



5. Illustrative examples and benchmarks 

We modeled a 2D freestanding slab irradiated from below on its 
sides and bottom by a Planck function. The slab is homogeneous 
and static with a vertical geometrical extension z ma x - 30 000 
km; its horizontal extension y max could take the respective val- 
ues: 100000, 30000, 10000, 5 000 and 1 000 km. Depth points 
are logarithmically spaced away from the boundary surfaces and 
the graphical representation we adopted compresses the central 
region and greatly expand the ar eas ne ar the boundaries. We 
have used the "set A" of Carlson (11963b with 3 points per oc- 
tant to describe the angular dependence of the radiation field and 
constant Doppler profiles. The temperature of the slab was fixed 
to 7=5 000 K and the gas pressure p g = 1 dyncm -2 . Finally, 
we adopted the standard benchmark models for multilevel atom 
problems proposed by Avrett ( 1968| see also Paletou & Leger 
120071) considering, in particular, its 3-level H i atomic model. 

The respective rates of convergence for the SOR and MG-2D 
multilevel iterative processes are displayed in Fig. [6] where we 
have plotted the maximum relative change on the level popula- 
tions (i.e., the oo-norm) from an iteration to another R c . The com- 
putation time for the MALI, GS, SOR and MG 2D-multilevel 
iterative processes are given in Tab. [T] for different grid refine- 
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— ID case 

- - y™*= 100 kkm 



0.8^ 




03 1 10 10 2 10 3 10 4 10 5 10 6 10 7 10 8 i. 5 io 9 10 8 10 7 10 6 10 5 10 4 10 3 10 2 10 1 

z [cm] 



Fig. 7. Vertical variations of the Ha source function (in units of 
B) along the symmetry axis of the 2D slab with z max = 30000 
km and different horizontal extensions y max ranging from 1 000 
km to oo (ID); the temperature of the atmosphere is 5000 K 
and the gas pressure p g = 1 dyncm -2 . The solid line represents 
S(Ha)/B variations computed in ID. Note that abscissae give 
geometrical positions computed downwards from the top sur- 
face up to mid-slab and then, symmetrically, upward from slab 
bottom. 



ments. We remind that a MG scheme is not only superior in it- 
eration numbers and computing time: it is also important to note 
that the convergence error C e which is defined by 

l\n(itr)-n(oo)\\ 
C e = max — (28) 

\ rt(oo) / 

is smaller than R c for MG whereas for methods such as MALI 
or SOR a small value of R c does not imply a small value of 
C e , which means that convergence is not necessarily achieved 
(Fabiani Bendicho et al. U997b . 

As shown in Fig. [7] where S(Ha) normalized to the external 
illumination is plotted as a function of the vertical line-center 
optical depth, the same variations as in ID (solid line) are re- 
covered along the vertical axis of symmetry of the 2D model 
which has the largest horizontal extension (i.e., 100000 km). 
For smaller geometrical slab widths and accordingly horizon- 
tal optical thickness, lateral radiative transfer effects take place 
and progressively affect the excitation within the slab. Note that 
for the smallest width (i.e., 1 000 km), we properly recover an 
almost constant value S/B - 0.5 consistent with optically thin 
conditions along the horizontal extension of the 2D slab. As first 
reported by Paletou (1997 ), we also recover here under which 
conditions 2D radiative transfer effects on the Ha source func- 
tion vertical variations can be significant; more generally, such 
effects are a priori expected for any other spectral lines of mod- 
erate optical thickness. 

Figures[8]are contour plots of the two excited levels of hydro- 
gen normalized to their LTE values obtained for a 30 000 km by 
5 000 km slab. They show both departures from LTE together 
with geometry effects within the slab atmosphere. However, 
since such data do not exist yet in the litterature while needs for 
multidimensional radiative modelling tools are more and more 
obvious, and in order to detail the information content of Fig. [8] 
about the populations distribution across the 2D slab, we found 
that Tab.|2]can also be highly valuable for benchmark purposes. 




Fig. 8. Contour diagram of (top) the second-level populations 
^2/^2 anc ^ (bottom) the third-level populations n^/n* in a grid 
with y max =5 kkm, z ma x =30 kkm, T=5 000 K and p g = 
1 dyncm -2 . The slab is illuminated from below with a Planck 
function and n* 2 and rf^ are LTE values. Horizontal and vertical 
axes are defined as described in Fig. [7] 

6. Conclusions 

We have given here details upon the implementation of GS/SOR 
iterative processes in 2D cartesian geometry, information which 
was unfortunately still missing in the astrophysical litterature. 
We also tested, for the first time, such 2D-GS/SOR iterative 
schemes with a two-level atom model against original analyti- 
cal results; a more comprehensive study, both in ID and in 2D, 
is being conducted and results will be published elsewhere. 

Concerning the modelling of illuminated freestanding slabs, 
even though we used here a quite simple atomic model, we found 
it to be a necessary stage not only to valid our numerical work 
but also to take the opportunity to deliver reliable 2D multilevel 
benchmark results; typical CPU usage numbers were also given, 
clearly in favour of the combination of SOR plus MG methods 
for complex radiative modelling. 

We anticipate that such numerical techniques and benchmark 
results will be of interest for the new radiative transfer codes 
currently in use or under development, not only for applications 
in solar physics but also for interstellar clouds (see e.g., Juvela & 
Padoan 2005), circumstellar environments with winds (see e.g., 
Georgiev et al. 2006) or accretion disks (see e.g., Korcakova & 
Kubat 2005) modelling for instance. 

Acknowledgements. Our warmest thanks go to Dr. Bernard Rutily for the origi- 
nal idea and fruitful discussions upon the analytical test presented here; we also 
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Table 2. Second-level (top) and third-level (bottom) populations for the Avrett (11968b Hi atomic model in a 2D grid of 163 points 
per direction with y max = 5000 km and z m ax = 30000 km together with 3 angles per octant and 8 frequencies ; the temperature of 
the atmosphere is 5000 K and the gas pressure p g = 1 dyn cm -2 . 
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296.091 
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295.636 


3 x 10 9 - 


10 


179.929 
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0.161275 


10 2 




2.39169 
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0.422406 
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2.59030 
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2.04905 


1.61663 


1.38949 
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2.65493 
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3.04375 
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2.33988 


2.66628 
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3.91497 


3.91163 
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2.66898 


2.58607 


2.34158 


2.66891 


3.53832 
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3.92572 


3.92885 


3.92452 
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10 8 
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3.56437 


3.92360 
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2.75824 


3.65778 
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1.5 x 10 9 


2.72068 


2.64248 


2.41682 


2.76842 


3.67139 


4.04146 


4.07319 


4.07500 


4.07460 


4.07610 


3 x 10 9 - 
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2.64845 


2.42479 


2.77894 


3.68546 


4.05695 


4.08873 
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4.07652 


3 x 10 9 - 
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2.86803 
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3 x 10 9 - 
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2.49389 
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4.19223 
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3 x 10 9 - 


10 4 
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2.89035 
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4.20084 
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4.22730 


4.23259 


4.23593 


3 x 10 9 - 


10 3 
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4.25128 


3 x 10 9 - 


10 2 


3.08898 
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3.45736 


4.22654 
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4.52946 


4.53266 


4.53348 


4.53832 


4.54137 


3 x 10 9 - 
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3.71314 


3.95438 


4.74064 


5.11335 


5.21097 


5.22648 


5.22770 


5.22826 


5.23210 


5.23452 


3 x 10 9 - 


- 1 


3.99111 


4.39634 


5.06123 


5.33812 


5.40867 


5.41987 


5.42076 


5.42127 


5.42483 


5.42707 



thank an anonymous referee for her/his valuable comments which helped us to 
clarify some technical points. 
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Appendix A: Test case for a 2D code using 1 D 
reference solutions 

We describe a test case for radiative transfer methods in 2D 
cartesian geometry with stationary media, using ID reference 
solutions, which are provided using an analytical method. For 
this purpose, the ARTY code is the numerical implementation, 
whose accuracy is better than 10" 10 , of exact analytical solu- 
tions, based on a mathematical method using the finite Laplace 
transform (Chevallier & Rutily 2005, Chevallier et al. 2003 and 
references therein). 

Our radiative model describes a 2D medium which can scat- 
ter in 3D and is infinite and homogeneous along the x-axis 
(-00 < x < +oo, < v < y max , < z < Zmax), thus quan- 
tities involved in the radiative transfer equation (RTE) do not 
depend on x. This medium is considered such that there is no 
incoming flux on its boundaries along the y- and z-axes. In order 
to compare this 2D case to ID solutions from ARTY, we con- 
sider here the 2D primary source to be an infinite line along the 
x-axis, located at the center of the slab, emitting isotropically, 
and the medium homogenous and isotropically scattering; the 
later is also monochromatic i.e., the RTE does not depend on the 
frequency (which will not be mentioned hereafter) as this is the 
case when we describe the continuum or a spectral line with the 
Milne profile, which is constant over any finite energy range and 
elsewhere. 

We write hereafter the RTE in 2D cartesian geometry, and 
we show how to compare this 2D solution integrated on the y- 
axis to a ID solution. Table I A. ll resumes some values of the ID 
solution at the surface z = 0. The RTE for our 2D model is 
(cf. Chandrasekhar 1950, Chap. I, Eq. (48) or Pomraning 1973, 
Eq. (2.60), without derivative over x though) 



dl dl 
sin#sin<£— (y,z, 0, (p) + cos0— (y,z,0,<p) = 
oy oz 



- x U(y,z,e,<p)-S(y,z)].. 



(A.l) 



where / is the specific intensity of the radiative field at (y,z) 
and in the direction (0, (p) of the unit vector n whose coordinates 
along x, y and z are sin cos (p, sin 6 sin (p and cos 0, respectively. 
X is the constant opacity of the homogeneous medium, and S is 
the unknown source function which can be written 



S(y,z) = S*(y,z) + viJ(y,z), 



(A.2) 



where S* describes the primary source function i.e., the direct 
known radiative field emitted by internal sources, tu - (1 - s) is 
the constant scattering coefficient of the homogeneous medium 
for simple scattering processes, usually called albedo, and / is 
the mean intensity of the radiative field defined as 



1 r»27T 

— \ d6 I t 

4?r Jo Jo 



dip sin 1(y, z, 0, tp) ; 



the primary source function is 

S (y,z) = -6(y- ——)5(z- — ), 
X 2 2 



(A.3) 



(A.4) 



where L is the luminosity per unit length along the x-axis. 
Dividing by x means that the source function is an emissivity 
divided by the opacity. In order to use ID solutions as a refer- 
ence, we must integrate the 2D solutions on y over [0, y max ] and 
on (p over [0, 2n]. We thus define new functions as 



/fc 0) = ±- f ymaX dy f d<pl(y, z, 0, ip) . 
*n Jo Jo 



(A.5) 



Similarly we define S(z) = S*(z) + mJ(z\ S*(z) = L/x$(z - 
Zmaxii), J(z) and the two successive moments, the radiative flux 
H(z) and the radiative pressure K(z) as 



[J,H,K](z) 



1 r 

2 Jo 



I(z, 0)[1, cos 0, cos 2 0] sin OdO . 



(A.6) 



Integrating over y and <p, and using the symmetry property valid 
for (p e [0, 7r] : I(y max , z, 0, <p) - 1(0, z,0,n + (p), due to the central 



primary source, Eq. (A.l } becomes 



sine r . 

sm cpl(y 

n Jo 

-X[I(z,6)-S(z)], 



dl 

, z, 0, (p)d(p + cos 0— (z, 6) = 
oz 



(A.l) 



where the integral is nul only for 6 = or n\ note that this simpli- 
fication is fictitious as, even for these angles, the source function 
depends on the mean intensity which depends on the boundaries 
due to the angular integration. This problem is not classical and 
we need to let y max —> +oo in order to suppress this term i.e., the 
radiation of the primary source is nul at the infinite and Eq. ( A.7 ) 
then reduces to the well-known ID equation: 



dl 

jd—(z,/d) = -xU(z,fi) - S (z)] , 
oz 



(A.8) 



where yu = cos 0. 

Equation (A.8 ) is usually expressed in optical depth coordi- 
nates t(z) = f w x(z')dz' = x(Zmax-z). We do not write the RTE, 
but the primary source function becomes S *(r) = L6(r - r max /2) 
due to the Dirac transformation S(z) = X^(xz)- Accordingly our 
2D primary source function becomes 



S*(Ty,T z ) =XLS(Ty ~ T -^)S(T Z ~ ^p), 



(A.9) 



where r y = x^max ~ i) and r z = x(z ma x ~ z). In order to simplify 
the test of a 2D code with a ID reference solution, the values 
L - 1 and^ = 1 shoul d be used. 

We give in Table |A.1| some values of the ID solution at 
the surface z = 0, for the source function, the specific in- 
tensity for the directions of the angular grid used in this pa- 
per, and its three first moments. When integrating all angles 
over the azimuthal angle (p, the 10-points per octant angular 
quadrature resume to a 4-points per quadrant, i.e. [/, H, K](z) = 
Ei = i,4 Wj[l,//i,/i?]/(z,/ii) for such a case where there is no 
incoming flux. The four directions \i[ are 0.95118969679, 
0.78679579496, 0.57735025883, 0.21821789443 and the in- 
tegration weights Wi are 0.063490696251, 0.091383516788, 
0.12676086649, 0.21836490929, respectively. It is interesting to 
note that, using the reference solutions, the angular quadrature 
for J, H and K will lead to a relative error equal to 0.8%, 0.3% 
and 0.4% respectively. 
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Table A.l. Reference solutions from the ARTY code at the sur- 
face z = for our test case with L = 1, x - 1> Zmax = 100 and 
s = 0.01 (see the text for the values of //,•). 



ARTY results 


5(0) 


2.710704655 x 10" 


-4 


J(0) 


2.738085511 x 10" 


-4 


H(0) 


1.600980711 x 10" 


-4 


K(0) 


1.130399095 x 10" 


-4 


/(0,//i) 


7.837047273 x 10" 


-4 




6.965481933 x 10" 


-4 


/(0,yU 3 ) 


5.880635905 x 10" 


-4 


/(0,// 4 ) 


4.026985767 x 10" 


-4 



